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We report a study of the process e + e“ —»• (D*D*)° 7r° using e + e _ collision data samples with 
integrated luminosities of 1092 pb~ 4 at yfs = 4.23 GeV and 826 plW 1 at y/s = 4.26 GeV collected 
with the BESIII detector at the BEPCII storage ring. We observe a new neutral structure near the 
(D*D*)° mass threshold in the 7r° recoil mass spectrum, which we denote as Z c ( 4025)°. Assuming a 
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Breit-Wigner line shape, its pole mass and pole width are determined to be (4025.5l 2 '7±3.1) MeV/c 2 
and (23.0 ± 6.0 ± 1.0) MeV, respectively. The Born cross sections of e + e“ —> Z c (4025)°n° —4 
(D*D*)° 7t° are measured to be (61.6 ± 8.2 ± 9.0) pb at yfs = 4.23 GeV and (43.4 ± 8.0 ± 5.4) pb at 
y/s = 4.26 GeV. The first uncertainties are statistical and the second are systematic. 


PACS numbers: 14.40.Rt, 13.25.Gv, 13.66.Be 

Recent discoveries of new charmoniumlike states that 
do not fit naturally with the predictions of the quark 
model have generated great experimental and theoret¬ 
ical interests [l]. Among these so-called “ XYZ ” par¬ 
ticles are charged states with decay modes that clear¬ 
ly demonstrate a structure consisting of at least four 
quarks, including a cc pair. The first charged charmoni¬ 
umlike state Z( 4430) + was discovered by Belle (2|. LHCb 
confirmed the existence of this state. Belle determined 
its spin-parity to be 1 + [1], which is supported by a 
new result from LHCb[||. Recently, the BESIII collab¬ 
oration observed four charged Z c states, ^' c (3885) ± @], 
2’ c (3900) ± Z c { 4020) ± [7], and Z c (4025) ± @] , pro¬ 
duced in e + e“ —4 n T Zf. The observed decay chan¬ 
nels are Z c (3900) ± -A t ^J/ip, Z c { 3885) ± -4 (DD*)*, 
Z c (4020) ± -4 7 ^hc, and Z £ (4025) ± -A (£>*£>*)*. These 
states are close to the DD* or D*D* threshold. The 
.Zc (3900) ± was also observed by Belle Q and with 
CLEO-c data E3- 

So far, the nature of these new states is still elu¬ 
sive. Interpretations in terms of tetra-quarks, molecules, 
hadro-charmonium, and cusp effects have been pro¬ 
posed Em Searching for their neutral partners in 
experiment is of great importance to understand their 
properties, especially to identify their isospin properties. 
Previously, based on CLEO-c data, evidence of a neutral 
state Z c (3900)° decaying to ir°J/ip [20j was reported. 
Recently, two neutral states, Z c (3900)° and Z c (4020)°, 
were discovered in their decays ZM3900) 0 —4 n^J/ip 
and Z c (4020)° -A n°h c by BESIII gj [22|. These can 
be interpreted as the isospin partners of the Z c (3900) ± 
and Z c (4020) ± . Analogously, it is natural to search for 
the neutral partner of the Z c (4025) ± @] in its decay to 
(. D*D*)°. 

In this Letter, we report a search for the neutral part¬ 
ner of the Z c (4025) ± through the reactions e + e _ — 4 
D*°D*°(D* + D*~)ir°, as the charged Z c (4025) ± [§| cou¬ 
ples to ( D*D*) ± and has a mass close to the ( D*D*) ± 
mass threshold. We denote the investigated final state 
products as (D*D*)°tt °, where D* refers to D*° or D* + , 
and D* stands for their antiparticles. A partial recon¬ 
struction method is applied to identify the (D*D*)°n° fi¬ 
nal states. This method requires detection of a D and a D 
originating from D* and D* decays of D* — 4 Dn and D'y, 
and the 7r° from the primary production (denoted as the 
bachelor 7r°). The data sample analyzed corresponds to 
e + e _ collisions with integrated luminosities of 1092 pb^ 1 
at y/s = 4.23 GeV and 826 pb -1 at Vs. = 4.26 GeV Q 
collected with the BESIII detector |24[ at the BEPCII 
storage ring @. 

BESIII is a cylindrically symmetric detector, which 


from inner to outer parts consists of the following com¬ 
ponents: a Helium-gas based multilayer drift chamber 
(MDC), a time-of-flight counter (TOF), a CsI(Tl) crystal 
electromagnetic calorimeter (EMC), a 1-Tesla supercon¬ 
ducting solenoid magnet and a 9-layer RPC-based muon 
chamber system. The momentum resolution for charged 
tracks in the MDC is 0.5% at a momentum of 1 GeV/c. 
The energy resolution for photons in EMC with energy 
of 1 GeV is 2.5% for the center region (barrel) and 5% 
for the rest of the detector (endcaps). For charged par¬ 
ticle identification (PID), probabilities C(h) for particle 
hypotheses h = tt or K are evaluated based on the nor¬ 
malized energy loss dE/dx in the MDC and the time of 
flight in the TOF. More details on the BESIII spectrom¬ 
eter can be found in Ref. [24} . 

To optimize data-selection criteria, understand back¬ 
grounds and estimate the detection efficiency, we simu¬ 
late the e + e~ annihilation processes with the KKMC algo¬ 
rithm [2fj , which takes into account continuum process¬ 
es, initial state radiation (ISR) return to ip and Y states, 
and inclusive production. The known decay rates 
are taken from the Particle Data Group (PDG) [27 ( and 
decays are modeled with evtgen [2Ej]. The remaining 
decays are simulated with the LUNDCHARM package [2£j] . 
The non-resonant, three-body phase space (PHSP) pro¬ 
cesses e + e~ —4 are simulated according to uni¬ 

form distributions in momentum phase space. We as¬ 
sume that .Z c (4025) 0 has spin-parity of 1 + by consid¬ 
ering the measurements of other Z resonances [1, [|| 
and the signal process e + e“ — 4 Z c ( 4025)°7r° followed by 
Z c (4025)° -A ( D*D*)° proceeds in S waves. The D* 
is required to decay inclusively according to its decay 
branching ratios from PDG [27|. The D + is required to 
decay into K~Tr + -ir + while D° is required to decay into 
K~ 7r + , K~ 7r + 7r° and K~tt + tt + tt~. These decay modes 
are the ones used to reconstruct D mesons [30]. All sim¬ 
ulated MC events are fed into a GEANT4-based [3lJ soft¬ 
ware package, taking into account detector geometry and 
response. 

The charged tracks of K~ and are reconstruct¬ 
ed in the MDC. For each charged track, the polar an¬ 
gle 9 defined with respect to the e + beam is required 
to satisfy |cos0| < 0.93. The closest approach to the 
e + e _ interaction point is required to be within ±10 cm 
along the beam direction and within 1 cm in the plane 
perpendicular to the beam direction. A track is iden¬ 
tified to be a K(tt) when the PID probabilities satisfy 
C(K) > C(tt) {£{K) < £(7r)), according to the informa¬ 
tion from dE/dx and TOF. 

The 7 r° candidates are reconstructed by combining 
pairs of photons reconstructed in the EMC that are not 
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associated with charged tracks. For each photon, the en¬ 
ergy deposition in the EMC barrel region is required to 
be greater than 25MeV, while in the end-cap region, it 
must be greater than 50 MeV, due to the different detec¬ 
tor resolution and probabilities of reconstructing a fake 
photon. To suppress electronics noise and energy de¬ 
posits unrelated to the event, the EMC cluster time is 
restricted to be within a 700 ns window near the event 
start time. The invariant mass of any pair of photons 
M( 77 ) is required to within (0.120, 0.145) GeV/c 2 and 
constrained to the nominal 7 r° mass. The kinematics of 
the two photons are updated according to the constraint 
fit. 

We consider all possible combinations of selected 
charged tracks and 7 r° to form D candidates. The charged 
tracks from a D decay candidate are required to origi¬ 
nate from a common vertex. The Xvf of the vertex fit is 
required to satisfy x'vf < 100. We constrain the recon¬ 
structed masses of the final state particles to the corre¬ 
sponding D nominal masses and require XkfC^) f° r the 
kinematic fit to be less than 15 for the final states of D 
decays including charged tracks only, and less than 20 for 
the final state including 7 r°. We select signal event candi¬ 
dates which consist of at least one pair of DD candidates 
that do not share particles in the final state. If there is 
more than one pair of DD candidates in an event, only 
the one with the minimum Xkf(-^) + Xkf(^) kept for 
further analysis. 

We reconstruct the bachelor 7r° from the remaining 
photon showers that are not assigned to the DD pair. To 
further reject backgrounds, each photon candidate origi¬ 
nating from the bachelor 7 r° is required not to form a 7r° 
candidate with any other photon in the event. A mass 
constraint of the two photons to the 7r° nominal mass is 
implemented and the corresponding fit quality is required 
to satisfy Xkf( 7I ’°) < 20. To reject background for the 
bachelor 7 r° from D* —> Dn° decays, we require the Dtt° 
invariant mass to be greater than 2.02 GeV/c 2 . 

To identify the decay products of the signal process 
e + e _ —> D*D*n °, we plot the recoil mass spectra of 
Dtt° ( RM(Dir °)), as shown in Fig. [U The peaks around 
2 GeV/c 2 correspond to the process e + e _ —> DD*n° with 
a missing D*. Besides these peaks, we see clear bumps 
around 2.15 GeV/c 2 in data. These bumps are consis¬ 
tent with the MC simulations of the D*D*tt° final state. 
The peak position roughly corresponds to the sum of 
the mass of D* and the mass of a 7r, since the 7 r orig¬ 
inating from D* is soft and is not used in the compu¬ 
tation of the recoil mass. The backgrounds beneath the 
bumps are mostly from ISR production of D*D* process. 
Other processes, such as e + e _ —> £)*£)** — D*D*n°, are 
expected to be absent according to simulation studies. 
This is understandable because the process _Dg(2400) —»■ 
D*tt° is forbidden due to the conservation of spin-parity. 
Hj'(2420)° (D 2 ( 2460 )°) is narrow, and the sum of the 
mass of Z?j'(2420)° (£>2(2460)°) and D* is much larger 
than 4.26 GeV. To extract the signals, we keep events 
within the two-dimensional oval regions in the distribu- 
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FIG. 1. Distributions of RM(Dn°) at y/s = 4.23 GeV (a) 
and yfs = 4.26 GeV (b). Points with error bars are data and 
the shaded histograms represent the inclusive backgrounds in 
MC simulations. The soild line and the dashed line are the 
-27(4025)° signal shape and the PHSP shape with arbitrary 
normalization, respectively. The third row gives the scatter 
plot of RM(D-k°) versus RM(Dn°) at yfs — 4.23 GeV (c) and 
yfs = 4.26 GeV (d) , where the solid ovals indicate the signal 
regions. 

tions of RM(Dn°) and RM(Dtt°) shown in Fig. [ljc,d). 
We choose the specific dimensions due to different resolu¬ 
tions at different momentum phase spaces at two energy 
points. They are determined according to MC simula¬ 
tion. 

The selected events are used to produce the recoil 
mass distribution of the bachelor 7 r° ( RM(tt °)), shown 
in Fig. [2j We observe enhancements in the RM(ir°) dis¬ 
tribution over the inclusive backgrounds for both data 
samples, which can not be explained by three-body non¬ 
resonant processes. We assume the presence of an S’-wave 
Breit-Wigner resonance structure (denoted as i? c (4025) 0 ) 
with a mass-dependent width, using the form given in 
Ref. [12: 


M 2 -m 2 -i-m{T 1 (M) + T 2 (M))/c 2 F ’ 

and r fc (M) = / fc -r-^~ (k = 1,2). 

Pk M 

Here, k = 1 and 2 denote the neutral chan¬ 

nel -2T C (4025)° —> £)*°£)* 0 and the charged channel 
Z c ( 4025)° —> D* + D *~, respectively. /*, is the ratio of 
the partial decay width for channel k. M is the recon¬ 
structed mass, m is the resonance mass and F is the reso¬ 
nance width. p k (qk) is the D*(tt°) momentum in the rest 
frame of the D*D* system (the initial e + e _ system) and 
pi is the momentum of D* in the Z c (4025)° rest frame 
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at M = to. We assume that Z c ( 4025)° decay rates to the 
neutral channel and the charged channel are equal, i.e., 
fk = 0.5, based on isospin symmetry. 

We perform a simultaneous unbinned maximum like¬ 
lihood fit to the spectra of RM(tt°) at y/s = 4.23 and 

4.26 GeV. The signal shapes are taken as convolutions of 
the efficiency-weighted Breit-Wigner functions with res¬ 
olution functions obtained from MC simulations. The 
detector resolutions are 4MeV at yfs = 4.23 GeV and 

4.5 MeV at y/s = 4.26 GeV. Backgrounds are mod¬ 
eled with kernel-estimated non-parametric shapes [33j 
based on the inclusive MC, and their magnitudes are 
fixed according to the simulations, since the inclusive 
MC samples well describe the background. The shape 
of the PHSP process is adopted from MC simulations. 
We combine the data at yfs = 4.23 GeV and yfs = 

4.26 GeV together, as shown in Fig. [2] The fit de¬ 
termines m and T to be (4031.7 ± 2.1)MeV/c 2 and 
(25.9 ± 8.8) MeV, respectively. The corresponding pole 

position m P oie(-Zc(4025) 0 ) — j ; r p° 1 <d + ( 4 "-.>) 1 is calculated 
to be 

m pole (Z c (4025)°) = (4025.5t 2 ;°)MeV/c 2 , 

r pole (Z c (4025)°) = (23.0 ± 6.0) MeV. 

The significance with systematic errors is estimated by 
comparing the likelihoods of the fits with and without 
the Z c ( 4025)° signal component included. The likeli¬ 
hood difference is 2AlnL = 45.3 and the difference of 
the number of free parameters is 4. When the systematic 
uncertainties are taken into account with the assumption 
of Gaussian distribution, the significance is estimated to 
be 5.9cr. 

The Born cross section a(e + e~ —> Z c ( 4025)°7r° —» 
(. d*°D*° + D* + D*~)7 t°) is calculated from the equation 

77-sig 

a “ CUiBiei + f 2 B2£2){l + S){ 1 + <5 vac ) ’ 

where C is the integrated luminosity, £\ (£ 2 ) is the detec¬ 
tion efficiency of the neutral (charged) channel, /1 (/ 2 ) 
is the ratio of the cross section of the neutral (charged) 
channel to the sum of the both channels, B\ (£> 2 ) is the 
product branching fraction of the neutral (charged) D* 
decays to the final states we detected. (1 + S) is the 
radiative correction factor and (1 + 5 V ac ) is the vacuum 
polarization factor. From the simultaneous fit, we obtain 

69.5 ± 9.2 signal events at yfs = 4.23 GeV and 46.1 ± 8.5 
signal events at y/s = 4.26 GeV. (1+5) is calculated to be 
0.744 at 7s = 4.23 GeV and 0.793 at ^/i = 4.26 GeV to 
the second order in QED [3] , where the input line shape 
of the cross section is assumed to be the same as for 
e + e _ —> (D*D*) + tt ~, as extracted directly from BESIII 
data. (1 + 5vac) is given as 1.054 following the formula 
in Ref. i35i|. The efficiency £\ (£ 2 ) is determined to be 
1.49% (3.87%) at y/s = 4.23 GeV and 1.84% (4.37%) at 
y/s = 4.26 GeV. Thus, the cross sections are measured to 



RMl+KGeV/c 2 ) 


FIG. 2. Fits to (a) A fit to background, PHSP 

and Zc(4025)° signal process for the combination of all da¬ 
ta (main plot), and the two collision energies separately (in¬ 
set plots), (b) Fits using only the inclusive background and 
PHSP. Points with error bars are data, solid line is the sum 
of fit functions, dotted line stands for the Z c (4025)° signals, 
filled area represents inclusive backgrounds, and dash-dotted 
line is the PHSP process. 


TABLE I. Summary of systematic uncertainties on the 
Z c (4025)° resonance parameters and cross sections 04230 at 
y/s = 4.23 GeV and 04260 at 4.26 GeV. “■ • • ” means the un¬ 
certainty is negligible. The total systematic uncertainty is 
taken as the root of the quadratic sum of the individual un¬ 
certainties. 


Source 

m(MeV/c 

2 ) r(MeV) (74230 (%) 

<74260 (%) 

Tracking 



5 

5 

Particle ID 



5 

5 

7T° reconstruction 



4 

4 

Photon veto 



4.2 

4.2 

Mass scale 

2.6 




Detector resolution 

0.2 

0.1 

0.3 

0.5 

Backgrounds 

0.6 

0.2 

5.6 

5.4 

Oval cut 

1.5 

1.0 

4.2 

2.0 

Fit range 


0.1 

0.3 

0.5 

D* D* 7T° line shape 



6.0 

3.0 

Luminosity 



1 

1 

£?i and B 2 



6.5 

5.3 

Isospin violation 


0.2 

0.3 

0.2 

Vacuum polarization 



0.5 

0.5 

Total 

3.1 

1.0 

14.6 

12.5 


be (61.6 ± 8.2) pb and (43.4 ± 8.0) pb at yfs = 4.23 and 

4.26 GeV, respectively. The contribution of the PHSP 
process is found to be negligible according to the fit. 

Sources of systematic uncertainties in the measure¬ 
ment of the Z c (4025)° resonance parameters and cross 
sections are listed in Table Q] Uncertainties of tracking 
and PID are each 1% per track f3f|. The uncertainty 
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of the 7T° reconstruction efficiency is 4% [37} . We study 
the photon veto by fitting the recoil mass of Dtt° with 
and without this veto in selecting the control sample of 
e + e _ —> (D*D*)°tt° in data. The efficiency-corrected 
signal yields are used to extract the cross section, and 
the corresponding change is taken into account as the 
systematic error introduced by this requirement. The 
systematic uncertainties are determined to be 4.2% for 
both data samples. The mass-scale uncertainty for the 
Z c ( 4025)° mass is estimated with the mass shift (com¬ 
parison between the PDG nominal values and fit values) 
of RM(Dn°) in the control sample e + e - —> DDtt 0 and 
of RM{D) in the control sample of e + e _ —>■ DD. To be 
conservative, the largest difference of the two mass shifts, 
2.6 MeV/c 1 2 3 , is assigned as the systematic uncertainty due 
to the mass scale. The systematic uncertainty from back¬ 
grounds is estimated by leaving free the magnitudes in 
the fit and making different choices in non-parametric 
kernel-estimation of the background events to account for 
the limited precision in MC simulation jdV|. We change 
the oval cut criteria and take the largest difference as the 
systematic uncertainty. Since the line shape will affect 
the efficiency and (1 + 5), to evaluate the systematic un¬ 
certainties with respect to the input D*D*ir a line shape, 
we change its shape based on uncertainties of the ob¬ 
served D* + D*°n~ cross section. Branching fractions B\ 
and #2 are used in calculating the cross sections and the 
uncertainties of the world average results are included as 
part of the systematic uncertainty. 

Other items in Table U have only minor effects on the 
precision of the results. We change the fitting ranges in 
the RM(ir°) spectrum and take the largest difference as 
the systematic uncertainty. The uncertainties due to de¬ 
tector resolution are accounted for by varying the widths 
of the smearing functions. The uncertainty of integrated 
luminosity is determined to be 1% by measuring large 
angle Bhabha events [7]. We vary the ratio fk from 
0.4 to 0.6 to take into account potential isospin viola¬ 
tion between the neutral and charged processes. The 
corresponding changes are assigned as systematic uncer¬ 
tainties. The systematic uncertainty of the vacuum po¬ 
larization factor is 0.5% [35} . 

In summary, using e + e _ annihilation data at yfs = 
4.23 and 4.26 GeV, we observe enhancements in the 
7 r° recoil mass spectrum in the process e + e _ —> 
D*°D*°(D* + D*~) 7T°. Assuming that the enhancement is 


due to a neutral charmoniumlike state decaying to D*D* 
and it has spin-parity of 1 + , the mass and width of its 
pole position are determined to be m P oie(^c(4025)°) = 
(4025.5t4;?±3.1)MeV/c 2 and P po i e (.2 c (4025) 0 ) = (23.0+ 
6.0 ± 1.0) MeV, respectively. The Born cross section 
a(e + e~ ->■ Z c (4025)°7r° (D*°D*° + D*+D*~) tt°) is mea¬ 

sured to be (61.6 ± 8.2 ± 9.0) pb at yfs = 4.23 GeV 
and (43.4 ± 8.0 ± 5.4) pb at y/s = 4.26 GeV. Hence, 

we estimate the ratio !. n*\+ 

cr(e+e —>Z C (4025) + n D*) + 7r ) 


to be compatible with unity at y/s = 4.26 GeV, which 
is expected from isospin symmetry. In addition, the 
-Z c (4025)° has mass and width very close to those of the 
.Z c (4025) ± , which couples to (0*0*)^ Q. Therefore, the 
observed Z c (4025)° state in this Letter is a good candi¬ 
date to be the isospin partner of Z c ( 4025) ± . 
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